Prostorno - vremenska predikcija Normalnog vazdušnog pritiska koristeći tehnike geostatistike i mašinskog učenja
Student: Petar Bursać; Profesor: prof. Dr Milan Kilibarda, Aleksandar Sekulić
28 jul 2020
Uvod
Cilj ovog rada je analiza i predikcija dnevnih karata Normalnog vazdšnog pritiska (Sea level pressure) na području teritorije R. Srbije, na osnovu merenja dostupnih sa 65 mernih stanica OGIMET mreže. Prostorni domen predikcije definisan je pravilnim gridom, čija prostorna dispozicija odgovara prostornom rasporedu mernih stanica. Vremenski domen merenja obuhvata period od 10 godina (2009-2018 godine), sa dostupnim merenjima na dnevnom nivou, izraženih kao srednje vrednosti dnevnog prikupljanja podataka sa senzora. U radu je korišćen programski jezik R sa dostupnim paketima, kojima je omogućena implementacija Prostorno-vremenske predikcije korišćenjem tehnika geostatistike i mašinskog učenja (spacetime, gstat, GSIF, randomForest, caret, ranger).
Metodologija
Metodologija podrazumeva jasno definisan niz postupaka u cilju izrade karata Normalnog vazdušnog pritiska sa prostornom rezolucijom od 1km. Metodologija obuhvata korišćenje sledećih pristupa u oceni i predikciji posmatrane promenljive:
- Oridnary kriging - Obični kriging
- Regression kriging - Regresioni kriging (Ocena trenda linearnom regresijom, kvanitifikacija i modelovanje prostorne korelacije na osnovu reziduala i predikcija vrednosti reziduala)
- Machine learning techique Random Forest - Korićenje tehnika Mašinskog ucenja, tj. Random Forest modela u cilju ocene trenda, a zatim predikcija posmatrane promenljive na osnovu definisanog modela i prediktora.
Kao ocena kvaliteta predikcije izvršen je postupak krosvalidacije:
- Metodom Leave-one station out (LOOCV) - sukcesivno isključivanje pojedinačnih stanica iz modela i ponavljanje procesa predikcije na lokaciji posmatrane stanice.
Sračunate su sledeće mere ocene kvaliteta predikcije krosvalidacijom:
- Srednja kvadratna greška - \(\mathbf{RMSE}\) (Root Mean Square Error)
- Koeficijent determinacije - \(\mathbf{R^2}\)
- Lin-ov koeficijent korelacije saglasnosti
Analiza i istraživanje seta podataka
Analizom dostupnog seta podataka na raspolaganju su sledeći podaci:
- 65 mernih stanica OGIMET mreže stanica, sa vrednostima koordinata tačaka u WGS84 referentnom koordinatnom sistemu.
- Vremenska serija podataka koja odgovara periodu od 01.01.2009 - 31.12.1018 godine (10 godina - 3652 dana). Kao jedinica mere, definisan je [dan].
- Vrednosti posmatrane promenljive SLP - Sea Level Pressure u jedinicama mere [mb], sa ukupnim brojem vremenskih trenutaka, na svim stanicama, bez vrednosti opažanja - 23527.
| country | n |
|---|---|
| Albania | 3 |
| Bosnia and Herzegovina | 3 |
| Bulgaria | 4 |
| Croatia | 5 |
| Hungary | 7 |
| Macedonia, The Republic of | 7 |
| Montenegro | 7 |
| Romania | 5 |
| Serbia | 24 |
Prostorna dispozicija mernih stanica vazdušnog pritiska
Za potrebe analize seta podataka i istraživanja (sagledavanja lokacija mernih stanica, njihove udaljenosti i prostornog rasporeda, visina na kojima se nalaze), veoma je bitno u početnoj fazi rada sagledati prostornu dispoziciju seta podataka, kao i postaviti početne hipoteze u vezi njihovih lokacija. Na sledećem sledećoj karti moguće je sagledati prostornu dispoziciju 65 mernih stanica OGIMET mreže stanica.
Prostorna dispozicija mernih stanica Normalnog vazdusnog pritiska data je i putem web kartografskog prikaza, koji je kreiran korišćenjem R programskog paketa “mapview” [Appelhans et al., 2018]
Vremenska serija podataka
Analizom seta podataka, kreiran je grafički prikaz vremenske serije podataka. Na gornjem delu grafičkog prikaza data je vremenska linija vrednosti Normalnog vazdušnog pritiska, po godinama, sa tim što su po \(\mathbf{x}\) osi date vrednosti vremenskih trenutaka izraženih po Julianskom kalendaru (godina 2455000 po Julianskom kalendaru odgovara 2009 godini po Gregorijanskom kalendaru), Na donjem delu grafika dat je boxplot grafički prikaz Normalnog vazdušnog pritiska izraženog po godinama.
Analizom grafičkih prikaza ne uočavaju se velike promenu u intenzitetu pritiska, ako se vrši poređenje po godinama. Zaključuje se da su vrednosti u potpunosti slučajnog karaktera, iz statističkog ugla posmatranja. Najmanja izmerena vrednost zabeležena je 2015. godine.
Na sledećem grafiku je moguće ostvariti uvid u vremensku liniju podataka po stanicama:
Analizom podataka i uvidom u grafički prikaz vremesnkih serija podataka po stanicama, jasno se uočava da na pojedinim stanicama postoje vremenski trenuci bez registrovanih vrednosti vazdušnog pritiska. Stanica 13619 je u potpunosti bez merenja za dati vremenski period. U kasnijem nastavku rada dati vremenski trenuci bez merenja, kao i stanica, su uklonjeni iz seta podataka.
Histogram
Na sledećem grafičkom prikazu dat je histogram frekvencija merenih vrednosti vazdušnog pritiska, kao i kriva normalne rapodele koja odgovara datom uzorku. Uočava se da se podaci u potpunosti pokoravaju normalnoj raspodeli, sa srednjom vrednošću i merama odstupanja vrednosti slučajne promenjive od njenog matematičkog očekivanja - standardnom devijacijom. Kriva raspodele je blago pomerena, i postoje uzorci u setu podataka čiji intenzitet je izvan matematčkog očekivanja. Kasnijom obradom podataka potrebno je ostvariti detaljniji uvid u set podataka kao i izvršiti statisitčko testiranje hipoteza da u skupu podataka nema statistički značajnih grubih grešaka. Glavni preduslov za korišćenje tehnika predikcije - kriginga je da podaci, tj. uzorak merenja pripada normalnoj raspodeli, što je i ispunjeno.
Hovmoller grafički prikaz
Hovmöller plot - grafički prikaz je dvodimenzionalni prostorno-vremenski način vizuelizacije. Prostorni domen je predstavljen na jednoj osi (projektovan ili osrednjen), dok je po drugoj osi definisano vreme. Ovakav način vizuelizacije je tradicionalno primenjivan u naukama i oblastima koje se bave praćenjem stanja atmosfere, mora i okeana, kako bi se na veoma efikasan način predstavila - vizuelizovala propagacije entiteta. [Wikle et al., 2019] Pomoću ovakvog načina vuzuelizacije moguće je sagledati korelacije u prostorno-vremenskom domenu, tj. npr. geografske širine ili dužine i vremena, vremenskih trenutaka opažanja posmatrane promenljive. Na ovaj način moguće je ostvariti i preliminarni uvid u postajenje trenda koji se opisuje i koordinatama tačaka.
Analzom Hovmöller plot-ova po koordinatnim osama definisanog referentnog koordinatanog sistema (WGS84), kao i vremena - vremenskih trenutaka opažanja, uočava se da ne postoji znatna prostorno-vremenska korelacija, koja bi opisivala korelisanost Normalnog vazdušnog pritiska opažanog tokom vremena, kao i lokacija stanica (njihove geografske dužine i širine).
Empirijska prostorna sredina (Empirical Spatial Mean)
Veoma koristan način analize prostorno-vrmenskih podatka je računanje empirijske sredine. Ako se uzme u obzir da su dostupna opažanja u definisanom prostornom i vremenskom domenu, empirijska prostorna sredina za svaku pojedinu lokaciju - mernu stanicu je sračunata agregacijom vrednosti tokom vremena. Empirijsku prostornu sredinu je moguće vizuelizovati, kao što je prikazano na sledećem scatter plot-u (i Howmoller plot je jedan od načina vizuelizacije). Cilj ove analize je pre svega kvantifikacija prostorno-vrmenske zavisnosti, korelisanosti, kao i ocene i vrednovanja potencijalnog trenda u podacima.
Na datom grafičkom prikazu vrednosti pritiska u odnosu na geografsku širinu i dužinu, koje su dobijene agregacijom tokom vremena, jasno se uočava da postoji zavisnost između geografske širine i vazdušnog pritiska. Veća vrednost koeficijenta determinacije ovo potvđuje, ali još uvek nedovoljno, čime se postavlja hipoteza za dalja istraživanja i statisitičko testiranje.
Emprijska vremenska sredina (Empirical Temporal Means)
Empirijska vremenska sredina predstavlja agregaciju vrednosti promenljive na lokacijama od interesa tokom vremena, sa tim što sad agregacija se vrši u vremenskom domenu koji je posmatran, za razliku od empirijske prostorne sredine gde su lokacije mernih stanica bile od interesa. Agregacija podataka je sračunata tokom godina, meseca i dana, a dostupna je na sledećem grafičkom prikazu.
Graficki prikaz koji sadrzi srednje vrednosti normalnog vazdusnog pritiska iskazanom po danu, mesecu i godini
Uvidom u vremesnke serije podataka, osrednjene tokom definisanih vremenskih intervala, kao i prostorno sumiranih, jasno se uočava pojava paterna i uticaja doba godine (seasonal pattern) i drugih vremenskih prilika i atmosferskih promenljivih na vazdušni pritisak. Uočava se periodična promena raspona vrednosti na par godina, ali godišnje srednje vrednosti su približno jednakih vrednosti.
Prosta linearna regresija (Simple Linear Regression)
Cilj proste linearne regresije je opisivanje i kvatifikovanje veze izmedju dve neprekidne slucajne velicine. Ako je:
- \(\mathbf{y_i}\) - vrednost zavisne promenljive i-te opservacije u uzorku
- \(\mathbf{x_i}\) - vrednost prediktora (nezavisne promenljive) i-te opservacije u uzorku
Dati uzorak neće pripadati tačno lineranoj pravoj, ali dovoljno je da bude blizu sa nekim manjim odstupanjem, pa da linearna veza bude dovoljno objašnjenje. Dakle, model je:
\[ \begin{aligned} \mathbf{y_i=b_0+b_1x_i+e_i} \end{aligned} \]
Parameter \(\mathbf{b_0}\) se naziva intercept - presek, a \(\mathbf{b_1}\) slope - nagib. Ovaj model se naziva prosta linearna regresija (Simple Linear Regression). Ako je u pitanju linearna zavisnost jedne zavisne sa više nezavisnih promenljivih radi se o višestrukoj linearnoj regresiji (Multiple Linear Regression). Iako se koriste nazivi zavisna i nezavisna promenljiva, vrednosti \(\mathbf{x_i}\) se ne smatraju slučajnim u modelu. Slučajan je niz grešaka \(\mathbf{e_i}\), pa odatle i niz \(\mathbf{y_i}\).
Ako je dat dvodimenzioni uzorak (x,y) sa grafika raspršenja (scatter plot) dobija se informacija o postojanju i obliku zavisnosti.
Scatter plot sa koeficijentima korelacije:
Kod linearne regresije parametri modela su odabrani tako da je suma kvadrata reziduala najmanja. Postavlja se pitanje i koliko je ta minimalna suma mala, jer ona govori o tome u kojoj meri je linerani model odgovarajući, koliko je disperzije početnog skupa objašnjeno uzimanjem u obzir linerane veze sa prediktorom. Jedna mera koja govori o tome koliko je smanjeno srednje kvadratno odstupanje je koeficijent determinacije \(\mathbf{R^2}\). Koficijent determinacije je izražen na sledeći način:
\[ \begin{aligned} \mathbf{R^2 = RSS/SSTO = 1 - SSE/SSTO} \end{aligned} \]
gde je:
- SSR - Regression Sum of Squares - suma kvadrata reziduala regresije
- SSTO - Total Sum of Squares - ukupna suma kvadrata reziduala
- SSE - Error Sum of Squares - suma grešaka kvadrata reziduala
Ako je \(\mathbf{R^2 = 1}\) to znači da sve tačke \(\mathbf{y_i}\) pripadaju regresionoj pravoj. Ako je \(\mathbf{R^2 = 0}\) onda je ocenjena linija horizontalna, pa prediktor ne objašnjava varijabilnost u zavisnoj promenljivoj ništa bolje od srednje vrednosti. Nije moguće na osnovu \(\mathbf{R^2}\) porediti modele koji imaju i one koji nemaju intercept. Naziv \(\mathbf{R^2}\) potiče od toga što je ovako definisana vrednost jednaka +/â’ korenu koeficijenta korelacije između prediktora i zavisne promenljive. Vrednost \(\mathbf{R^2}\) se može često zloupotrebiti i pogrešno shvatiti, što je opisano sledećim primerima:
- \(\mathbf{R^2}\) meri linearnu vezu između datih promenljivih. Ako se dobije da je \(\mathbf{R^2 = 0}\) to znači da nema linearne zavisnosti, ali možda postoji neka nelinearna zavisnost.
- Velika vrednost \(\mathbf{R^2}\) ne mora da znači da je linearna funkcija najbolja veza. Ovo se može videti sa grafika.
- Jedna neobična vrednost u uzorku može puno da utiče na vrednost \(\mathbf{R^2}\).
- Velika vrednost \(\mathbf{R^2}\) ne znači prediktivnu moć.
Grid - domen predikcije
Za potrebe kreiranja predikcije vrednosti Normalnog vazdušnog pritiska nad područjem od interesa i u vremenskom domenu, koji je definisan kao dan, kreiran je grid - skup tačaka na pravilnom i jednakom rastojanju. Prostorna rezolucija grid-a je 1km, dok je vremenska rezolucija 1 dan, posmatrano za vremenski interval 01.01.2009-31.12.2018. Ciljni grid je kreiran kao objekat STF klase podataka, paketa “spacetime” (Spatio-Temporal Feature).
Prostorno-vremenski obični kriging (Spatio-temporal Ordinary Kriging)
U okviru ovog poglavlja dat je opis, metodologija i rezultati predikcije Normalnog vazdušnog pritiska kao dnevnih karata na području od interesa, obuhvaćenog mernim stanicama. Metodologija se zasniva na korišćenju Kriging tehnika predikcije u prostorno-vremenskom domenu.
Sva merenja - opažanja na površi Zemlje, iako je ponekad ignorisano, imaju prostorno - vremensku komponentu (referencu). Prostorno vremenska komponenta opažnja je određena sa minimum četri parametra:
- Geografska - prostorna lokacija (geografskim koordinatama ili u definisanom koordinatnom sistemu u projekciji).
- Visina na samoj površi Zemlje u odnosu na definisanu referentnu nivosku površ (visina).
- Vreme merenja - opažnja (godina, mesec, dan, sat, itd.).
- Prostorno - vremensku podršku - određenost (veličinu bloka pridruženih merenjima, vremenski interval merenja;).
Na osnovu prethodno definisanih komponenti moguće je definisati čitavi niz metodologije obrade i analize podataka u mnogim oblastima, a zajedničko za tu grupu nauka je da se analize definšu kao prostorno - vremenske analize podataka (Spatio - temporal Data Analysis). Grubo rečeno prostorno vremenske analize podataka se mogu posmatrati kao kombinacija dve osnovne nauke: geoinformatika i prostorno vremenska statistika. Druga reč koja se može upotrebiti za termin prostorne statistike, a isto tako dodeliti i vremensku dimenziju, je nauka koja se zove Geostatistika (Geostatistics).
Geostatistika je grana matematičke statistike koja se bavi analizom i interpretacijom geografski referenciranih podataka, tj. prostorno određenih i kojima je moguće dodeliti vremensku dimenziju. Geostatistika pre svega ima fokus na prostonim podacima koji su kontinualni (prostorna kontinualnost) i pritom se postavljaju osnovna pitanja na koja je potrebno da odgovori: kako se promenljiva ponaša u prostorno - vremenskom koordinatnom sistemu, šta je uzrok i kontrola varijacije, gde je potrebno utorkovati određenu pojavu da bi se opisala prostorna promenljivost, koliko uzoraka je potrebno, kolika je vrednost promenljive na nekoj novoj lokaciji, kolika je nesigurnost ocenjenih vrednosti.
Postavljena pitanja opisiju glavni zadatak na koji Geostatistika kao nauka je potrebno da odgovori. Sama metodologija je jasno definisana i osnovne hipoteze koje se postavljaju pred svaku promenljivu koja egzistira u prostoru su sledeće:
- Određivanje modela (Model estimation) - zaključci o parametrima modela.
- Predikcija - interpolacija vrednosti promenljive (Prediction) - zaključivanje o prediktovanoj vrednosti posmatrane promenljive.
- Testiranje hipoteza (Hypothesis testing).
Kao sinonim za Geostatističku interpolaciju može se reći da je to Kriging, kao statističku prostornu predikciju. Geostatističke metode su prvi razradili francuski geomatematičar Georges Matheron i južnoafricki rudarski inženjer D.G. Krige. Zato se ove metode u literaturi nazivaju kriging. U geodeziju je geostatistiku uveo austrijski geodeta Moritz (pod nazivom kolokacija) za potrebe odredivanja anomalija sile zemljine teže i otklona vertikala.
U oviru ovog rada date su osnovne teroijske postavke na kojima se zaniva Kriging kao metoda predikcije prostornih promenljivih. Dat je osvrt na kvantifikaciju prostorne strukture promenljivih, definisanje variograma i komponenti variograma, kao i koncepti Kriging predikcije. Eksperimentalni deo rada je urađen dostupnim programskim paketima, kojima je moguće izvršiti geostatističku analizu i interpretaciju u okviru programskog jezika R. Osnovne funkcionalnosti su dostupne putem paketa “gstat”, koji podržava i rad sa prostorno vremenskim promenljivima, čija implementacija je moguća kroz definisane klase paketa “spacetime”:
- STF - Spatio-Temporal Full grid layout - Prostorno - vremenski potpuni grid
- STS - Spatio-Temporal Sparse grid layout - Prostorno - vremenski redak grid
- STI - Spatio-Temporal Irregular layout - Prostorno - vremenski nepravilan grid
- STT - Spatio-Temporal Trajectory - Prostorno - vremenska trajektorija
Svaka od pomenutih klasa omogućuje i proširenje na -DF - Data Frame u cilju čuvanja osobina entiteta, atributa. Forma u kojoj je podaci mogu čuvati je:
- long formats - potpuna kombinacija prostornog i vremenskog domena
- space-wide - gde različite kolone predstavljaju različite lokacije u prostoru
- time-wide - gde različite kolone različite trenutke vremena
Evidentno da definisani pristup čuvanju i obradi podataka putem dostupnih paketa, u poptunosti predstavlja kompatibilno i efikasno rešenje za potrebe prostorno-vremenske predikcije koristeći Geostatističke metode.
Variogram - modelovanje i kvantifikacija prostorne korelisanosti
Interpolacija primenom geostatističkih metoda se izvodi u dva koraka:
- Kvantifikacija prostorne strukture površi koja se modelira (na osnovu ulaznih podataka), i
- Predikcije, tj. ocene vrednosti funkcije površi u zadatim tačkama.
Kvantifikacija prostorne strukture površi koja se modelira izvodi se empirijskim određivanjem kovarijacione funkcije, odnosno poluvariograma kojima se opisuje prostorna zavisnost vrednosti funkcije u tačkama površi. Osnovni problem kod praktične primene ovih metoda je upravo izbor odgovarajuće kovarijacione funkcije (funkcije poluvariograma), tj. izbor odgovarajućeg modela funkcije i empirijsko određivanje parametara te funkcije. Varijacije vrednosti promenljive koja se posmatra, najbolje se mogu opisati stohastičkim površima, gde se posmatrana pojava posmatra kao regionalizovana promenljiva. Vrednost promenljive \(\mathbf{Z}\) na lokaciji \(\mathbf{x}\) se tada posmatra kao realizacija jedne slučajne promenljive, \(\mathbf{Z(x)}\) je slučajna funkcija, a cela površ slučajni proces. Teorija regionalizovane promenljive pretpostavlja da se prostorne varijacije mogu razložiti na tri nkomponente:
\[ \begin{aligned} \mathbf{Z(x) = m(x) + ε′(x) + ε′′(x)} \end{aligned} \]
- determinističke varijacije (trend) - \(\mathbf{m(x)}\)
- slučajne, ali prostorno autokorelisane varijacije - \(\mathbf{ε′(x)}\)
- prostorno nekorelisani šum (očekivana vrednost šuma je \(\mathbf{0}\), a disperzija je \(\mathbf{σ^2)}\) - \(\mathbf{ε′′(x)}\)
Za slučajne promenljive u prostoru se kaže da su regionalizovane. Za ekstrakciju prostornih osobina regionalizovanih promenljivih možemo koristiti geostatističke mere. Jednom kvantifikovane, osobine regionalizovanih promenljivih mogu se koristiti u mnogo aplikacija. Geostatistička analiza koristi prostorne autokorelacione informacije u procesu kriging interpolacije. Fenomen da su geografski bliže pojave međusobno više korelisane nego one koje su dalje. Autokorelacija je statistička veza između merenih tačaka, čija korelacija zavisi od udaljenja i pravca pojedinačnih lokacija. Mi znamo iz posmatranja realnog sveta da prostorna autokorelacija postoji zato što generalno zapažamo da pojave koje su međusobno bliže su mnogo više slične nego one koje su udaljenije. Kako se udaljenje povećava, prostorna autokorelacija opada.
Variografija (Variography) je proces kvantifikacije modela prostorne zavisnosti koja proizilazi iz podataka i prostorne strukture. Za određivanje veličine interpolacije na specifičnoj lokaciji, kriging koristi postavljen model u variografiji, prostornu konfiguraciju podataka i vrednosti merenog uzorka tačaka oko lokacije za predikciju. Jedna od najvažnijih mera koja se koristi za razumevanje prostorne strukture regionalizovanih promenljivih je variogram, koji se može korsititi za odnos varijanse unutar dela prostora (i autokorelacije) između uzorka. Poluvarijansa daje nepristrasan opis razmere i obrasca/šablona/uzorka prostorne promenljive unutar regiona. Geostatistički model zahteva kvantitativno merenje prostorne korelacije (spatial corelation) radi ocene i simulacija i naviše korišćen alat kojim se meri, kvantifikuje prostorna korelacija je poluvariogram. Variogram je grafikon kojim se rastojanja prikazuju u formi korelacije. Na osonovu seta ulaznih podataka, merenja iz realnog sveta moguće je dobiti eksperimentalni variogram. Osnovne karakteristike su da na kraćim rastojanjima korelacija je veća, dok povećanjem rastojanja korelacija opada. Na određenim rastojanjima među podacima ne postoji korelacija.
Dve bitne pretpostvke je potrebno uzeti u obzir:
- γ - poluvarijansa je funkcija rastojanja - pretpostvka stacionarnosti (Stationarity assumption).
- Funkcija samo rastojanja, ali ne i pravca - pretpostvka izotropije (Isotropy assumption).
U slučaju da pravac ulazi u obzir, da prostorne razlike u određenim pravcima naglašene u odnosu na druge pravce, dolazi do pojave anizotropije (anistropy assumption). Ako postoji način moguće je odkloniti trendom. U slučaju pojave anizotropije definišu se azimuti i zatim definišu direkcioni variogram ili kompletni rasterski variogram (zonalni).
Definisana terminologija i pretpostavke su u potpunosti proširive u domenu vremena, tako da sada u slučaju Prosotorno-vremenskih podataka, potrebno je i vreme uzeti u obzir. Definicija poluvarijanse se dopunjuje u pogledu vrmemna, tako da sada poluvariogram dobija treću dimenziju kojom se opisuje i međusobna zavisnost podataka u vremenu, a ne samo u prostoru. U narednom delu biće detaljnije objašnjen koncept prostorno-vremenske zavisnosti, kao i način kvantifikacije i izrade empirijskog prostorno-vremenskog poluvariograma na osnovu merenja vazdušnog pritiska.
Transformacija podataka u koordinatni sistem u projekciji
Koordinate lokacija mernih stanica su izražene u elipsoidnim koordinatama. U cilju izrade variograma i kvantifikacije prostorne korelisanosti, potrebno je da rastojanja budu izražena u metrima, a ne u stepenima. U cilju toga izvršena je transformacija koordinata mernih stanica iz geografskog-elipsoidnog koordinatnog sistema WGS84 u koordinatni sistem u projekciji (UTM - Universal Mercator Projection - 34N zona), koji je dat sledećom notacijom u obliku EPSG koda (EPSG:32634):
\[ \begin{aligned} \mathbf{+init=epsg:32634 +proj=utm +zone=34 +datum=WGS84} \end{aligned} \] \[ \begin{aligned} \mathbf{+units=m +nodefs +ellps=WGS84 +towgs84=0,0,0} \end{aligned} \]
Variogram
Na sledećem grafičkom prikazu dat je prostorno-vremenski variogram opažanja vazdušnog pritiska:
Empirijski prostorno-vremenski variogram
Nakon kreiranja prostorno-vremenskog variograma, koji opisuje set podataka i zavisnost - korelaciju između merenja na različitim lokacijama i različitim vremenskim trenucima, potrebno je izvršiti fitovanje modela - određivanje matematičke funcije i njenih parametara koji na najbolji mogući način opisuju set podataka. Fitovanje variograma se postiže određivanjem inicijalnih parametara funkcije na osnovu vizuelnog pregleda podataka. Konačna empirijska funkcija, kao i njeni parametri se određuju automatski kroz više iteracija. R programski paket “gstat” nudi više modela za fitovanje variograma prostorno-vremenskih podataka: separable, product sum, metric, sum metric, and simple sum metric. Više informacija o svakom modelu ponaosob se može pronaću u radu [Graler et al., 2015]. Za modelovanje prostorno-vremenskog variograma vrednosti merenja vazdušnog pritiska, korišćn je sum metric model, koji podrazumeva uključivanje prostorne i vremenske komponente putem modela kovarijansi - poluvariograma, kao i zajedničke komponente sa vrednovanjam parametra anizotropije. Prostorno vremenska funkcija kovarijanse data je na sledeći način: \[ \begin{aligned} \mathbf{Csm (h, u) = Cs (h) + Ct (u) + Cjoint (sqrt(h^2 + (k*u)^2))} \end{aligned} \]
gde je:
- \(\mathbf{h}\) - rastojanje u prostornom domenu [m]
- \(\mathbf{u}\) - rastojanje u vremenskom domenu [dan]
- \(\mathbf{C}\) - matrica kovarijansi (s - prostorne, t - vremenske i joint - prostorno-vremenske komponente)
- \(\mathbf{k}\) - parametar anizotropije
U cilju određivanja inicijalnih vrednosti parametara variograma potrebno je zadati početne vrednosti parametra psill - parcijalni prag i range - domet. Vrednost parcijalnog praga je približno jednaka vrednosti varijanse seta podataka. Za početnu vrednost dometa može uzeti 1/4 vrednosti dijagonale minimalnog obuhvatnog pravougaonika područja od interesa obuhvaćenog stanicama, sa tim da se ta vrednost uzima samo za prostornu komponentu.
Srednja vrednost varijanse vrednosti Normalnog vazdušnog pritiska na osnovu vrednosti po stanicama [mb^2]
| Variance_s | units |
|---|---|
| 50.5537 | mb^2 |
Srednja vrednost varijanse vrednosti Normalnog vazdušnog pritiska na svim stanicama za posmatrani vremenski period [mb^2]
| Variance_t | units |
|---|---|
| 52.384 | mb^2 |
Srednja vrednost varijanse vrednosti Normalnog vazdušnog pritiska u prostorno-vremenskom domenu [mb^2]
| Variance_j | units |
|---|---|
| 51.4689 | mb^2 |
Prikaz empirijskog prostorno-vremenskog variograma, fitovanog korišćenjem sum-metric modela. Model je dobijen korišćenjem funkcija vgmST() i fit.StVariogram(), R programskog paketa “gstat”. Približna vrednost koeficijenta anizotropije je sračunata empirijski korišćenjem funkcije estiStAni(). Za fitovanje modela variograma korišćena je sferna funkcija kod sve tri komponente: prostorne, vremenske i zajedničke - prostorno-vremenske.
Mean Squared Error (MSE) - srednja kvadratna greška optimizacije modela variograma po metodi najmanjih kvadrata
| Mean Squared Error | units |
|---|---|
| 3.419 | mb^2 |
| model | psill | range |
|---|---|---|
| Nug | 0.0000 | 0.0 |
| Sph | 1.4154 | 188193.3 |
| model | psill | range |
|---|---|---|
| Nug | 5.3006 | 0.000 |
| Sph | 35.3674 | 8.019 |
| model | psill | range |
|---|---|---|
| Nug | 0.000 | 0.0 |
| Sph | 1.683 | 188193.3 |
| anisotropy | units |
|---|---|
| 5.052 | km/day |
Predikcija
Predikcija nad setom podataka u cilju dobijanja dnevnih karata Normalnog vazdušnog pritiska je urađena korišćenjem Običnog prostorno-vremenskog kriginga. Kriging tehnike su implementirane u okviru R programskog paketa “gstat”, uključijući i rad sa podacima u formatu podataka definisanim paketom “spacetime”. Prostorni domen predikcije je definisan gridom sa rastojanjem tačaka od 1km, dok je vremenski domen predikcije definisan vremenskim periodom dostupnog seta podataka, za period od 2009-2018 godine.
Krosvalidacija
Ocenu kvaliteta modela i predikcije kontinualne promenljive je moguće oceniti unutrašnjim i spoljašnjim merama kvaliteta. Kao unutrašnje mere kvaliteta, koriste se koeficijent determinacije lineranog modela i varijansa kriging predikcije. Kod spoljašnjih mera ocena kvaliteta modela predikcije, može se koristiti nezavisan set podataka kojim se model ocenjuje tokom ili nakon izgradnje modela. Kada nije dostupan set podataka za validaciju modela, moguće je iskoristiti isti set podataka, ali isključivanjem pojedinih prostornih i/ili vremesnkih entiteta - merenja iz izgradnje modela. Takva metoda ocene kvaliteta i validacije modela predikcije naziva se \(\mathbf{krosvalidacija}\).
Krosvalidacijom se vrši ocena modela predikcije na taj način što se za izgradnju modela koristi potpun set podataka, a zatim prilikom predikcije u zadatom prostorno-vremenskom domenu jedan i/ili više merenja na jednoj ili više lokacija uzoraka se ne koriste. Koriste se sledeće metode krosvalidacije:
- \(\mathbf{LOOCV}\) - Leave-one-out cross validation - isključivanje jedne stanice ili entiteta u vremenskom domenu tokom predikcije.
- \(\mathbf{5-FOLD}\) - Krosvalidacija sa isključivanjem 5 merenja u prostornom i/ili vremenskom domenu tokom predikcije ili deljenjem seta podataka na 5 grupa - foldova od kojih se jedan koristi za validaciju seta podataka, dok se preostali koriste za izgradnju modela.
- \(\mathbf{5-FOLD-neasted-CV}\) - Ugnježdena krosvalidacija sa 5 foldova, gde se podrazumeva kreiranje 5 grupa podataka, gde prilikom svake iteracije se jedna od grupa koristi za validaciju a preostale za izgradnju modela, i tako sukcesivno u što više nivoa i kreiranja podgrupa (sub-folds).
Date metode krosvalidacije nisu primenljive za modele koji ne podrazumevaju teoriju geostatistike, tj. teoriju prostorne korelacije. U ovom radu je korišćena prva metoda krosvalidacije - LOOCV u prostornom domenu sa pojedinačnim isključivanjem stanica iz predikcije. Korišćen je identičan model - prostorno-vrmenski variogram dobijen na osnovu potpunog skupa merenja. Sa obzirom da su dostupne 64 merne stanice sa opažanjima od 10 godina - ukupno 3652 dana, dobijeno je ukupno 233728 (64x3652) karata predikcije. Na osnovu dobijenih podataka sračunate su mere kvaliteta krosvalidacije: ukupna srednje kvadratna greška predikcije, koeficijent determinacije i Lin-ov koeficijent korelacije skladnosti.
Srednja kvadratna greška - \(\mathbf{RMSE}\) (Root Mean Square Error) predikcije, dobijena krosvalidacijom nad setom podataka, metodom LOOCV:
| RMSE | units |
|---|---|
| 1.0283 | mb |
Koeficijent determinacije \(\mathbf{R^2}\) predikcije, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| R2 | units |
|---|---|
| 0.9798 | / |
Koeficijent korelacije saglasnosti objedinjuje mere preciznosti i tačnosti da bi se utvrdilo koliko daleki posmatrani podaci odstupaju od linije savršene skladnosti (tj. pravca x=y na grafiku rasipanja posmatranih promenljivih - merene i preditovane vrednosti). Lin-ov koeficijent se povećava smanjenjem rastojanja između glave ose podataka do linije savrešene saglasnosti (tačnost podataka) i tačnosti podataka o njegovoj smanjenoj glavnoj osi (preciznost podataka). Lin-ov koeficijent korelacije saglasnosti, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| rho | units |
|---|---|
| 0.9898 | / |
Analizom rezultata krosvalidacije u vremenskom domenu, nakon agregacije u prostornom domenu, moguće je ostvariti uvid u ocenu kvaliteta kao i opseg vrednosti po mesecima. Očigledno je da postoji zakonitost u podacima, kojom se opisuje manji opseg vrednosti, kao i greška predikcije u letnjim mesecima. U zimskim mesecima postoji veća varijabilnost seta podataka, kao i da je greška predikcije veća.
| MESEC | RMSE [mb] | MAX [mb] | MIN [mb] | MEAN [mb] | RANGE [mb] |
|---|---|---|---|---|---|
| Jan | 1.28 | 1041.4 | 971.3 | 1018.84 | 70.1 |
| Feb | 1.28 | 1041.5 | 973.8 | 1016.64 | 67.7 |
| Mar | 1.43 | 1061.6 | 967.4 | 1016.28 | 94.2 |
| Apr | 1.24 | 1061.8 | 989.6 | 1014.91 | 72.2 |
| Maj | 0.91 | 1028.3 | 994.0 | 1013.99 | 34.3 |
| Jun | 0.82 | 1044.5 | 995.2 | 1013.79 | 49.3 |
| Jul | 0.71 | 1026.3 | 998.4 | 1013.64 | 27.9 |
| Avg | 0.68 | 1026.4 | 1003.0 | 1015.14 | 23.4 |
| Sep | 1.00 | 1077.8 | 998.3 | 1016.64 | 79.5 |
| Okt | 0.79 | 1040.2 | 991.6 | 1019.17 | 48.6 |
| Nov | 0.90 | 1043.7 | 987.8 | 1018.90 | 55.9 |
| Dec | 0.98 | 1041.8 | 992.5 | 1021.41 | 49.3 |
Agregacijom seta podataka u prostornom domenu, moguće je ostvariti uvid u grešku predikcije po stanicama, na kojima je opažan vazdušni pritisak za dati vremenski period.
Prediktori (Covariates, Explanatory variables)
Digitalni model terena (DEM - Digital Elevation Model)
Kao jedan od prediktora u cilju predikcije vrednosti Normalnog vazdušnog pritiska, korišćen je digitalni model terena na području od interesa. Korišćen je DEM nastao kao produkt SRTM satelistske misije (Shuttle Radar Topography Mission). Vazdušni pritisak opada sa visinom i to približno jedan milibar za svakih 10 metara povećanja visine. Razlog opadanju vazdušnog pritiska sa visinom je jednostavan. Prvo, penjanjem na veću visinu skraćuje sa vazdušni stub, pa onaj deo vazdušnog stuba koji ostaje ispod ne vrši pritisak. Drugo, sa visinom vazduh postaje sve ređi, a time i lakši. Već na 5.000 metara visine pritisak vazduha iznosi oko polovinu vrednosti pritiska na morskom nivou. Na 16.000 metara pritisak je već deset puta manji nego na morskom nivou, a na 30.000 – 100 puta manji. Zbog opadanja vazdušnog pritiska sa visinom meteorološke stanice, na različitim nadmorskim visinama, stanice će izmeriti u istom trenutku različite vrednosti pritiska. Radi toga vazdušni pritisak se određuje za morski nivo što omogučava upoređivanje i analizu podataka o vazdušnom pritisku.
Temperatura (vremenska serija dnevnih karata temperature vazduha 2009-2018)
U cilju predikcije Normalnog vazdušnog pritiska, kao prediktori kod korišćenja Prostorno - vremenskog Regresionog kriginga i Random Forest metode mašinskog učenja, u cilju ocene trenda posmatrane promenljive, iskorišćene su karrte temperature (srenjih, minimalnih i maksimalnih vrednosti) za posmatrani domen predikcije. Merene vrednosti su dostupne za svih 65 mernih OGIMET stanica za posmatrani vremenski period od 2009-2018 godine. Na osnovu merenih podataka dostupnih stanica, karte temmperatura sa prostornom rezolucijom od 1km i vremenskom rezolucijom od 1 dana, su dobijene korišćenjem Prostorno - vremenskog regresionog kriginga. Vrednost trenda je ocenjena na osnovu globalnog (STRK_global) modela [Kilibarda et al., 2014], gde su kao prediktori korišćeni geometrijski temperaturni trend, DEM i TWI podaci na području od interesa. Detaljna analiza i ocena predikcije temperatura je data u dokumentu TEMP_CV_RMarkdown.html.
Vremenska serija opаžanja temperatura na stnicama
Karte predikcije srednjih vrednosti temperatura - TMEAN
Karte predikcije mmaksimalnih vrednosti temperatura - TMAX
Karte predikcije minimalnih vrednosti temperatura - TMIN
Prostorno-vremenski Regresioni kriging (Spatio-temporal Regression Kriging)
U okviru ovog poglavlja dat je opis, metodologija i rezultati predikcije Normalnog vazdušnog pritiska primenom Prostorno - vremenskog Regresionog kriginga. Domen predikcije u prostornom smislu je područje od interesa obuhvaćeno mernim stanicama sa prostornom rezolucijom od 1km. Vremenski domen predikcije je period od 10 godina (2009 - 2018 godine) sa vremesnkom rezolucijom od 1 dana. Metodologija se zasniva na korišćenju tehnike Regresionog Kriginga u prostorno-vremenskom domenu. Vrednost trenda je ocenjena primenom Višestruke linearne regresije (MLR - Multiple Linear Regression), gde su kao prediktori korišćene dnevne karte srednjih, maksimalnih i minimalnih vrednosti temperatura, kao i digitalni model terena. Koeficijenti trenda, kao i vrednosti reziduala su ocenjeni na osnovu MLR modela, a zatiom je kreiran variogram na osnovu ocenjenih reziduala. Izvršena je ocena i fitovanje empirijskog prostorno-vremenskog variograma (sum-metric variogram, putem R paketa ‘gstat’).
Ocena trenda - linearni regresioni model
Model kreiran Višestrukom linearnom regresijom, gde su kao prediktori korišćene dnevne karte srednjih, maksimalnih i minimalnih vrednosti temperatura, kao i digitalni model terena.
##
## Call:
## lm(formula = slp ~ tmean + tmax + tmin + dem, data = df_slpTemp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.301 -3.462 0.401 3.840 65.298
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.016e+03 4.102e-02 24778.85 <2e-16 ***
## tmean -5.758e-01 1.282e-02 -44.92 <2e-16 ***
## tmax 6.146e-01 6.956e-03 88.36 <2e-16 ***
## tmin -4.553e-01 8.463e-03 -53.80 <2e-16 ***
## dem -2.661e-03 8.410e-05 -31.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.508 on 212805 degrees of freedom
## Multiple R-squared: 0.1928, Adjusted R-squared: 0.1927
## F-statistic: 1.27e+04 on 4 and 212805 DF, p-value: < 2.2e-16
Variogram - modelovanje i kvantifikacija prostorne korelisanosti
Na sledećem grafičkom prikazu dat je prostorno-vremenski variogram reziduala Normalnog vazdušnog pritiska na osnovu trenda ocenjenog linearnim modelom i definisanim prediktorima:
Empirijski prostorno-vremenski variogram
Nakon kreiranja prostorno-vremenskog variograma, koji opisuje set podataka i zavisnost - korelaciju između reziduala na različitim lokacijama i različitim vremenskim trenucima, potrebno je izvršiti fitovanje modela - određivanje matematičke funcije i njenih parametara koji na najbolji mogući način opisuju set podataka. Fitovanje variograma se postiže određivanjem inicijalnih parametara funkcije na osnovu vizuelnog pregleda podataka. Konačna empirijska funkcija, kao i njeni parametri se određuju automatski kroz više iteracija. R programski paket “gstat” nudi više modela za fitovanje variograma prostorno-vremenskih podataka: separable, product sum, metric, sum metric, and simple sum metric. Više informacija o svakom modelu ponaosob se može pronaću u radu [Graler et al., 2015]. Za modelovanje prostorno-vremenskog variograma vrednosti reziduala normalnog vazdušnog pritiska na osnovu ocenjenog trenda linearnim modelom, korišćn je sum metric model, koji podrazumeva uključivanje prostorne i vremenske komponente putem modela kovarijansi - poluvariograma, kao i zajedničke komponente sa vrednovanjam parametra anizotropije.
U cilju određivanja inicijalnih vrednosti parametara variograma potrebno je zadati početne vrednosti parametra psill - parcijalni prag i range - domet. Vrednost parcijalnog praga je približno jednaka vrednosti varijanse seta podataka. Za početnu vrednost dometa može uzeti 1/4 vrednosti dijagonale minimalnog obuhvatnog pravougaonika područja od interesa obuhvaćenog stanicama, sa tim da se ta vrednost uzima samo za prostornu komponentu.
Srednja vrednost varijanse vrednosti reziduala Normalnog vazdušnog pritiska na osnovu vrednosti po stanicama [mb^2]
| Variance_s | units |
|---|---|
| 65.7747 | mb^2 |
Srednja vrednost varijanse vrednosti reziduala Normalnog vazdušnog pritiska na svim stanicama za posmatrani vremenski period [mb^2]
| Variance_t | units |
|---|---|
| 67.2002 | mb^2 |
Srednja vrednost varijanse vrednosti reziduala Normalnog vazdušnog pritiska u prostorno-vremenskom domenu [mb^2]
| Variance_j | units |
|---|---|
| 66.4874 | mb^2 |
Prikaz empirijskog prostorno-vremenskog variograma, fitovanog korišćenjem sum-metric modela. Model je dobijen korišćenjem funkcija vgmST() i fit.StVariogram(), R programskog paketa “gstat”. Približna vrednost koeficijenta anizotropije je sračunata empirijski korišćenjem funkcije estiStAni(). Za fitovanje modela variograma korišćena je sferna funkcija kod sve tri komponente: prostorne, vremenske i zajedničke - prostorno-vremenske.
Mean Squared Error (MSE) - srednja kvadratna greška optimizacije modela variograma po metodi najmanjih kvadrata
| Mean Squared Error | units |
|---|---|
| 3.3154 | mb^2 |
| model | psill | range |
|---|---|---|
| Nug | 0.0000 | 0.0 |
| Sph | 2.6674 | 139217.1 |
| model | psill | range |
|---|---|---|
| Nug | 3.4738 | 0.0000 |
| Sph | 38.8430 | 7.8004 |
| model | psill | range |
|---|---|---|
| Nug | 0.0070 | 0.0 |
| Sph | 4.4299 | 351774.5 |
| anisotropy | units |
|---|---|
| 111263.6 | km/day |
Predikcija
Predikcija nad setom podataka u cilju dobijanja dnevnih karata Normalnog vazdušnog pritiska je urađena korišćenjem Regresionog prostorno-vremenskog kriginga. Kriging tehnike su implementirane u okviru R programskog paketa “gstat”, uključijući i rad sa podacima u formatu podataka definisanim paketom “spacetime”. Prostorni domen predikcije je definisan gridom sa rastojanjem tačaka od 1km na području od interesa, dok je vremenski domen predikcije definisan vremenskim periodom dostupnog seta podataka, za period od 2009-2018 godine.
Krosvalidacija
U ovom radu je korišćena metoda krosvalidacije - LOOCV u prostornom domenu sa pojedinačnim isključivanjem stanica iz predikcije. Korišćen je identičan model - prostorno-vremenski variogram dobijen na osnovu potpunog skupa podataka - reziduala na osnovu trenda ocenjenog višestrukom linearnom regresijom. Sa obzirom da su dostupne 65 merne stanice sa opažanjima od 10 godina - ukupno 3652 dana, dobijeno je ukupno 233728 (64x3652) karata predikcije. Na osnovu dobijenih podataka sračunate su mere kvaliteta krosvalidacije: ukupna srednje kvadratna greška predikcije, koeficijent determinacije i Lin-ov koeficijent korelacije skladnosti.
Srednja kvadratna greška - \(\mathbf{RMSE}\) (Root Mean Square Error) predikcije, dobijena krosvalidacijom nad setom podataka, metodom LOOCV:
| RMSE | units |
|---|---|
| 1.8257 | mb |
Koeficijent determinacije \(\mathbf{R^2}\) predikcije, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| R2 | units |
|---|---|
| 0.938 | / |
Lin-ov koeficijent korelacije saglasnosti, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| rho | units |
|---|---|
| 0.9695 | / |
Analizom rezultata krosvalidacije u vremenskom domenu, nakon agregacije u prostornom domenu, moguće je ostvariti uvid u ocenu kvaliteta kao i opseg vrednosti po mesecima. Očigledno je da postoji zakonitost u podacima, kojom se opisuje manji opseg vrednosti, kao i greška predikcije u letnjim mesecima. U zimskim mesecima postoji veća varijabilnost seta podataka, kao i da je greška predikcije veća.
| MESEC | RMSE [mb] | MAX [mb] | MIN [mb] | MEAN [mb] | RANGE [mb] |
|---|---|---|---|---|---|
| Jan | 1.89 | 1041.4 | 971.3 | 1018.84 | 70.1 |
| Feb | 1.93 | 1041.5 | 973.8 | 1016.64 | 67.7 |
| Mar | 1.90 | 1037.3 | 967.4 | 1016.63 | 69.9 |
| Apr | 1.98 | 1061.8 | 989.6 | 1014.82 | 72.2 |
| Maj | 1.93 | 1028.3 | 994.0 | 1014.00 | 34.3 |
| Jun | 1.76 | 1044.5 | 995.2 | 1013.99 | 49.3 |
| Jul | 1.71 | 1026.3 | 998.4 | 1013.81 | 27.9 |
| Avg | 1.78 | 1026.4 | 1003.0 | 1015.16 | 23.4 |
| Sep | 1.70 | 1077.8 | 998.3 | 1016.33 | 79.5 |
| Okt | 1.79 | 1040.2 | 991.6 | 1019.08 | 48.6 |
| Nov | 1.77 | 1043.7 | 987.8 | 1018.60 | 55.9 |
| Dec | 1.75 | 1041.8 | 992.5 | 1021.37 | 49.3 |
Agregacijom seta podataka u prostornom domenu, moguće je ostvariti uvid u grešku predikcije po stanicama, na kojima je opažan vazdušni pritisak za dati vremenski period.
Prostorno-vremenska predikcija koristeći Random Forest tehniku Mašinskog učenja (Spatio-temporal Random Forest)
U okviru ovog poglavlja dat je opis, metodologija i rezultati predikcije Normalnog vazdušnog pritiska primenom tehnika mašinskog učenja - Random Forest modela. Domen predikcije u prostornom smislu je područje od interesa obuhvaćeno mernim stanicama sa prostornom rezolucijom od 1km. Vremenski domen predikcije je period od 10 godina (2009 - 2018 godine) sa vremesnkom rezolucijom od 1 dana. Metodologija se zasniva na korišćenju Random Forest modela u prostorno-vremenskom domenu, u cilju predikcije posmatrane promenljive korišćenjem dodatnih nezavisnih promenljivih - prediktora u cilju occene trenda. Pored prediktora kojima se opisuje - vrši ocena trenda Normalnog vazdušnog pritiska (dnevne karte srednjih, maksimalnih i minimalnih temperatura i Digitalnog modela terena), korišćeni su i dodatni prediktori u cilju uključenja prostorne i vremenske komponente. Primenjena je metodologija koja se zasniva na korišćnju prostornih prediktora - Euklidskih rastojanja do lokacija merenja - mernih stanica. Metodologija je data i opisana u radu [Hengl et al., 2018]. Vremenska komponeta (na osnovu predložene metodologije) je uključena sa dva dodattna preidiktora: “cdate” - kumulativni dan posmatrane vremenske serije i “doy” - (“day of the year”) dan u godini. Ocena tačnosti modela je izvršena računanjem ukupne srednje kvadratne greške predikcije, koeficijenta determinacije i Lin-ovog koeficijenta korelacije skladnosti, nakon postupka krosvalidacije metodom LOOCV - stanice.
U cilju poređenja i ocene kvaliteta predikcije datom metodologijom [Hengl et al., 2018], izvršeno je testiranje modela sa i bez uključivanja prediktora kojima se opisuje vremenska komponenta na predloženi način. Predikcija u prostorno vremenskom domenu posmatrne promenljive primenom tehnnike mašinsog učenja - Random Forest se zasniva na sledećim koracima:
- Izbor prediktora i preklop posmatrane promenljive sa istima.
- Ocena trenda Random Forest modelom - treninng modela, sa izbor vrednosti parametara kojima se definiše model.
- Ocena znnačajnosti preiktora - variable importance.
- Krosvalidacija sukcesivnim isključivanjem pojedinačnih stanica u cilju ocene modela.
- Predikcija posmatrane u definisanom prostorno-vremenskom domenu na osnovu izbora finalnog modela.
Definisana metodologija je implementirana u programskom jeziku R putem paketa ‘ranger’, ‘caret’, ‘GSIF’, ‘spacetime’, ‘raster’, ‘ggplot2’, kao i ‘tidyr’ familije paketa.
Dodatni prediktori
U cilju uključivanja prostorne i vremnske komponente u model predikcije Random Forest metodom, izvršeno je kreiranje dodatnog seta podataka.
Prediktori - prostorna komponenta
Prediktori kojima je opisana prostorna komponenta, uzet u obzir prostorni položaj stanica i rastojanje do svake pojedinačno su “buffer-distance layers” [Hengl et al., 2018]. Kreirani lejeri predstavlju rasterske lejere gde svaki piksel sadrži vrednost rastojanja do posmatrane stanica. Na ovaj način uključena su Euklidska rastojanja - efekat geografske blizine u proces predikcije (slično kao kod Kriging tehnika predikcije).
Prediktori - vremesnka komponenta
Preiktori kojima je opisana vremenska komponenta su “cdate” i “doy”. Prediktor “cdate” predstavlja kumulativni dan za datu vremensku seriju podataka, dok prediktor “doy” (“day of the year” - dan u godini) predstavlja definisan redni broj dana u okviru posmatrane vremenske serije. Datim prediktorima se uključuje “rastojanje” u vremsnkom domenu. Promenljiva “doy” je bitna radi predstavljanja efekta sezonalnosti, dok je promenljiva “doy” - kumulativni dan bitna u cilju uključivanja dugoročnog trenda.
slp_sp@data$cdate = floor(unclass(as.POSIXct(as.POSIXct(paste(slp_sp@data$time), format=“%Y-%m-%d”)))/86400)
slp_sp@data$doy = as.integer(strftime(as.POSIXct(paste(slp_sp@data$time), format=“%Y-%m-%d”), format = “%j”))
Ocena trenda - variable importance
U cilju predikcije Normalnog vazdušnog pritiska izvršeno je kreiranje Random Forest modela na osnovu dostupnih prediktora u posmatranom prostorno-vremenskom domenu. Pristup je zasnovana na metodologiji dostupnoj u radu [Hengl et al., 2018]. Razlog korišćenja Random Forest tehnike mašinskog učenja u cilju predikcije posmatranne kontinualne promenljive, kakav je Normalni vazdušni pritisak, se ogleda u sledećem: putem Random Forest tehnike moguće je u model uključiti ne-linearne veze između zavisne i nezavisnih promenljivih i drugo, moguće je detekotavati i oceniti značajne interakcije između nezavisnih promenljivih. Sa obzirom da je Normalni vazdušni pritisak izrazito kontinualna promenljiva, gde se u prostornom domenu dešavajju veoma male promene na bliskim rrastojanjima, ne-linearne veze sa nezavisnih predikotorima nisu došle do izražaja u oceni trenda i poboljšanju modela. Random Forest tehnika je posebno atraktivna za izgradnju multivarijabilnih modela za prostornu predikciju, koji se mogu koristiti kao “baze znanja” u geonaukama. Neki od nedostataka Random Forest modela za prostornu predikciju su eksponencijalno povećanje intenziteta proračuna sa povećanjem trening podataka i prediktora, kao i visoka osetljivost kvaliteta predikcije na kvalitet ulaznih podataka. Ključ uspeha Random Forest tehnike se ogleda u kvalitetu trening podataka za izgradnju modela. Isto tako kvalitet prostornog uzorkovanja - lokacija merenih vrednosti posmatrane promenljive (radi minimiziranja problema s ekstrapolacijom i bilo kakvom vrstom pristranosti u podacima) i kvaliteta validacije modela (kako bi se osiguralo da se tačnost ne postiže prekomernim uklapanjem) posebno su važni [Hengl et al., 2018].
Dodatno poređenje izvršeno je uključivanjem dodatnih prediktora (“cdate”, “doy”) kojima je uključena vremenska komponenta.
Trening modela
Trening modela nad dostupnim setom podataka izvršen je korićnjem programskog jezika R i paketa “ranger”.
U cilju ocene modela, potrebno je izvršiti izbor parametara kojima se definiše izgradnja Random Forest modela. Kao jedan od načina izbora parrametara, moguće je izvršiti optimizaciju izbora parametara, tako da je skup parametara kojima se dobija minimalna greška modela, najoptimalniji izbor. Definisani su sledeći parametri:
- mtry: broj promenljivih koje se mogu podeliti na svakom čvoru. Važno je da vrednost mtry parametra bude veća od p/2, gde je p broj nezavisnih promenljivih (dok je optimalna vrednost 80% od ukupnog broja p).
- num.trees (number of trees): broj stabala.
- sample.fraction: podela uzoraka merenja na svakom stablu (vrednost > 0.6 se koristi kod regresije).
- min.node.size: minimalna veličina terminalnog čvora.
Model bez prediktora “cdate” i “doy”
Ranger result
Call: ranger(formula = formulaSLP, data = slp_sp@data, mtry = 55, min.node.size = 2, sample.fraction = 0.95, num.trees = 150, seed = 1, quantreg = TRUE, importance = “impurity”)
Type: Regression
Number of trees: 150
Sample size: 212810
Number of independent variables: 69
Mtry: 55
Target node size: 2
Variable importance mode: impurity
Splitrule: variance
OOB prediction error (MSE): 35.7317
R squared (OOB): 0.3190538
Model sa prediktorima “cdate” i “doy”
Ranger result
Call: ranger(formula = formulaSLP, data = slp_sp@data, mtry = 57, min.node.size = 2, sample.fraction = 0.95, num.trees = 150, seed = 1, quantreg = TRUE, importance = “impurity”)
Type: Regression
Number of trees: 150
Sample size: 212810
Number of independent variables: 71
Mtry: 57
Target node size: 2
Variable importance mode: impurity
Splitrule: variance
OOB prediction error (MSE): 1.118662
R squared (OOB): 0.9786814
Variable importance
Kao jedan od rezultata treninga modela je vrednost “variable importance” (značajnost promenljive). Ovim pristupom se direktno meri značajnost promenljive u izgradnji modela, posmatrajući kako nasumično premeštanje (random re-shuffling) (čuvajući distribuciju promenljive) svakog prediktora utiče na performanse modela.
Model bez prediktora “cdate” i “doy”
Model sa prediktorima “cdate” i “doy”
Krosvalidacija
Model bez prediktora “cdate” i “doy”
Srednja kvadratna greška - \(\mathbf{RMSE}\) (Root Mean Square Error) predikcije, dobijena krosvalidacijom nad setom podataka, metodom LOOCV:
| RMSE | units |
|---|---|
| 5.7473 | mb |
Koeficijent determinacije \(\mathbf{R^2}\) predikcije, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| R2 | units |
|---|---|
| 0.3705 | / |
Lin-ov koeficijent korelacije saglasnosti, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| rho | units |
|---|---|
| 0.5278 | / |
Analizom rezultata krosvalidacije u vremenskom domenu, nakon agregacije, moguće je ostvariti uvid u ocenu kvaliteta kao i opseg vrednosti po mesecima.
| MESEC | RMSE [mb] | MAX [mb] | MIN [mb] | MEAN [mb] | RANGE [mb] |
|---|---|---|---|---|---|
| Jan | 5.81 | 1040.9 | 981.7 | 1016.19 | 59.2 |
| Feb | 5.57 | 1040.7 | 989.3 | 1016.06 | 51.4 |
| Mar | 5.41 | 1041.2 | 979.2 | 1015.89 | 62.0 |
| Apr | 5.88 | 1041.4 | 979.2 | 1015.96 | 62.2 |
| Maj | 5.75 | 1042.3 | 978.8 | 1017.33 | 63.5 |
| Jun | 6.21 | 1061.8 | 988.4 | 1015.90 | 73.4 |
| Jul | 4.94 | 1044.5 | 971.3 | 1016.90 | 73.2 |
| Avg | 6.56 | 1041.5 | 988.8 | 1019.22 | 52.7 |
| Sep | 4.74 | 1043.7 | 986.7 | 1016.92 | 57.0 |
| Okt | 6.21 | 1077.8 | 986.0 | 1017.36 | 91.8 |
| Nov | 5.49 | 1040.8 | 991.6 | 1015.98 | 49.2 |
| Dec | 6.09 | 1041.0 | 967.4 | 1015.59 | 73.6 |
Model sa prediktorima “cdate” i “doy”
Srednja kvadratna greška - \(\mathbf{RMSE}\) (Root Mean Square Error) predikcije, dobijena krosvalidacijom nad setom podataka, metodom LOOCV:
| RMSE | units |
|---|---|
| 1.2856 | mb |
Koeficijent determinacije \(\mathbf{R^2}\) predikcije, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| R2 | units |
|---|---|
| 0.9685 | / |
Lin-ov koeficijent korelacije saglasnosti, dobijen krosvalidacijom nad setom podataka, metodom LOOCV:
| rho | units |
|---|---|
| 0.9836 | / |
Analizom rezultata krosvalidacije u vremenskom domenu, nakon agregacije, moguće je ostvariti uvid u ocenu kvaliteta kao i opseg vrednosti po mesecima.
| MESEC | RMSE [mb] | MAX [mb] | MIN [mb] | MEAN [mb] | RANGE [mb] |
|---|---|---|---|---|---|
| Jan | 1.11 | 1040.7 | 985.7 | 1016.43 | 55.0 |
| Feb | 1.92 | 1041.4 | 984.0 | 1015.87 | 57.4 |
| Mar | 1.27 | 1041.0 | 978.8 | 1016.37 | 62.2 |
| Apr | 1.05 | 1040.7 | 979.2 | 1015.46 | 61.5 |
| Maj | 1.42 | 1041.8 | 973.8 | 1016.40 | 68.0 |
| Jun | 1.51 | 1040.5 | 967.4 | 1016.47 | 73.1 |
| Jul | 1.23 | 1077.8 | 988.8 | 1017.14 | 89.0 |
| Avg | 1.23 | 1061.8 | 989.8 | 1018.09 | 72.0 |
| Sep | 1.03 | 1043.7 | 986.7 | 1018.44 | 57.0 |
| Okt | 1.07 | 1040.7 | 971.3 | 1016.27 | 69.4 |
| Nov | 1.19 | 1040.2 | 975.3 | 1016.36 | 64.9 |
| Dec | 1.20 | 1061.6 | 986.8 | 1016.05 | 74.8 |
Predikcija - finalni model
U okviru ovog poglavlja date su karte predikcije Normalnnog vazdušnog pritiska dobijene predikcijom primennom tehnika Mašinskog učenja. Predikcija je izvršena primenom finalnog modela kojim se opisuje zavvisnost posmatrane i nezavisnih promenljivih - prediktora.
Analizom grafika Variable importance - značajnosti promenljivih, uočava se da model sa najboljim performansama (najboljom ocenom kvaliteta nakon postupka krosvalidacije), koji se dobija kada su u model uključeni prediktori kojima se opisuje vremenska komponenta (cdate i doy). Dati prediktori su od velike značajnosti u preidikciji Normalnog vazdušnog pritiska, a pre svega opisivanju vremenske zavisnosti i trenda koji je prisutan u vremenskom domenu (sezonalnost). Prediktor “cdate” omogućava kreiranje RF modela kojim je moguće vrednovati različite prostorne paterne za svaki dan ponaosob, sa obzirom i da su merenja na stanicama različita od dana do dana. Dodatnom analizom grafika značajnosti promenjivih, uočava se da veliku ulogu u opisu trenda Normalnog vazdušnog pritiska imaju prediktori temprature za posmatrani prostorno-vremenski domen. Očigledno je da postoji visoka korelacija i zavisnost između Normalnog vazdušnog pritiska i temperatura, sa obzirom da su obe promenljive kontinualne i u značanoj vezi iz ugla meteorološko-klimatoloških parametara. Prediktor kojim je dat Digitalni model terena, dolazi do izražaja u modelu bez eksplicitno date vremenske kommponente, ali njihovim uključivanjem ovaj prediktor, koji je konstantan u vremenskom domenu, gubi na značajnosti i interakciji u cilju predikcije vazdušnog pritiska.
Ocena kvaliteta modela i predikcije je data putem sračunatih parametara nakon postupka krosvalidacije. Tačnost predikcije dostiže vrednost koeficijenta determinacije \(\mathbf{R^2 = 0.37}\) u modelu bez i \(\mathbf{R^2 = 0.97}\) u modelu sa prediktorima kojim je uključena vremenska komponenta. Vrednost srednje kvadratne greške je u najboljem slučaju \(\mathbf{RMSE = 1.29mb}\), što ukazuje na visoku prediktivnu moć i tačnost predikcije Normalnog vazdušnog pritiska u prostorno-vremenskom domenu, Očigledno je da značajnu interakciju i cilju predikcije posmatrane promenljive imaju i “buffer-distance layers” preiktori, putem kojih je uključena prostorna komponenta u model (prostorni položaj i rastojanja između mernih stanica). Finalni model u cilju predikcije Normalnog vazdučnnnog pritiska se sastoji iz prediktora putem kojih je uključena prostorna i vremenska komponenta, kao i Digitalni model terena i srednje, masksimalne i minimalne vrednosti temperature za posmatrani prostorno-vremenski period. Prostorni domen predikcije je definisan gridom sa rastojanjem tačaka od 1km na području od interesa, dok je vremenski domen predikcije definisan za period od 2009-2018 godine sa vremenskom rezolucijom od 1 dana.
Analizom karata predikcije (dat je primer za tri dana u okviru posmatranog vremenskog perioda), uočava se značajnost prediktora “buffer-distance layers” koji imaju na predikciju u prostornom domenu. Uočavaju se koncetrični krugovi oko pojedinih lokacija stanica, a iz mogućeg razloga značajne interakcije pojedinog prediktora i/ili nedostataka stanica i merenja na pojedinim područjima. Normalni vazdušni pritisak je izrazito kontinualna promenljiva, ali isto tako uočava se i trend Digitalnog modela terena u prostornom domenu.
Zaključak
U okviru ovog rada prikazane su tri tehnike prostorno vremenske predikcije za potrebe predikcije Normalnog vazdušnog pritiska: Obični kriging, Regresioni kriging i Random Forest tehnika mašinskog učenja. U cilju ocene kvaliteta predikcije i poređenja performansi, sračunate su definisane mere kvaliteta nakon postupka Leave-one-out krosvalidacije. Analizom rezultata dolazi se do zaključka da metoda Običnog kriginga daje najbolju ocenu predikcije Normalnog vazdušnog pritiska, koji je izrazito kontinualna promenljiva.
| Metoda | R2 | RMSE | rho |
|---|---|---|---|
| Obični kriging | 0.98 | 1.03 | 0.99 |
| Regresioni kriging | 0.94 | 1.83 | 0.97 |
| Random Forest | 0.97 | 1.29 | 0.98 |
Analizom narednog grafičkog priloga, vremenske serije opažanaih vrednosti i prediktovanih za izabranu meteorološku stanicu, uočava se da su prediktovane vrednosti veoma bliske opažanim. Prediktovane vrednosti putem tehnike Običnog kriginga su najpribližnije merenim, dok prediktovane vrednosti putem tehnike Regresionog kriginga imaju najveće rasipanje vrednosti.
Metoda Random Forest daje bolje rezultate u odnosu na Regresioni kriging, a verovatnno zahvaljući dodatnim prediktorima koji su ušli u model (“cdate” i “doy”). Putem tih dodatnih prediktora je data vremenska komponenta, dok je kod regresionog kriginga sračunat prostorno-vremenski variogram i na taj način pored prostorne, izvršena kvantifikacija interakcije zavisne i nezavisnih promenljivih i u vremenskom domenu. Prednosti Random Forest metode u odnosu na Kriging tehnike ogledaju se u tome da nisu potrebne striktne statističke pretpostavke o raspodeli i stacionarnosti posmatrane promenljive, model je moguće proširiti prediktorima različitih vrsta i osobina. Očigledno je da prediktori uzeti u obzir, kao što su karte srednjih, maksimalnih i minimalnih temperatura i Digitalni model terena, veoma verno opisuju trend Normalnog vazdušnog pritiska, što se može zaključiti i iz rezultata predikcije.
Normalni vazdušni pritisak je meteorološki - klimatološki element, čije je poznavanje od velike važnosti, kako za oblasti geonauka tako i za klimatološko fenomene koji su od velike važnosti za život ljudi na planeti Zemlji. Pokazano je da se primenom geostatističkih i tehnika mašinskog učenja mogu veoma kvaliteno dobiti karte posmatrane promenljive u visokoj prostornoj i vremenskoj rezoluciji. Pokazano je da primena programskog jezika R i dostupnih paketa za rad sa prostorno - vremenskim podacima umnogome ima velikih prednosti u analizi podataka, automatizaciji procesa, vizuelizaciji kao i u dostupnim matematičkim i statističkim algoritmima, metodama i modelima.